function M = compute_data_moments_mod(X_monthly,Nu)
 
% Periods 1-6: pre shock
% Periods 7-14: post shock

% First equation: adoption rate
X                                                       = X_monthly(7:end,:);
TimeZeroX                                               = repmat(X_monthly(1,:),[8,1]);
dX                                                      = X - TimeZeroX;
TimeShock                                               = ones(size(dX));
TimeShock(1:3,:)                                        = 0;
%InitialCondition                                        = X_monthly(6,:)-X_monthly(1,:);
InitialCondition                                        = X_monthly(1,:);
InitialCondition                                        = repmat(InitialCondition,[size(dX,1) 1]);
InitialCondition                                        = (InitialCondition - mean(InitialCondition(:)))/std(InitialCondition(:));
Z1                                                      = [ones(size(dX(:))) TimeShock(:) InitialCondition(:) TimeShock(:).*InitialCondition(:)];
Y1                                                      = dX(:);
Mb1                                                     = Z1\Y1;
resids_sq                                               = (Y1 - Z1*Mb1).^2;

% Second equation (residuals from second equation)
Z2                                                      = [ones(size(dX(:)))];
Y2                                                      = resids_sq;

% Third equation: cross-district adoption variance
VarXT                                                   = var(X_monthly(7:end,:)-TimeZeroX,[],2);
IndicPost                                               = ones(size(VarXT));
IndicPost(1:3)                                          = 0;
%Z3                                                      = [ones(size(VarXT)) IndicPost];
Z3                                                      = [IndicPost ones(size(VarXT))];
Y3                                                      = VarXT;

% Fourth equation: within-district adoption variance
VarX                                                    = var(X_monthly(7:end,:)-repmat(X_monthly(1,:),[8,1]),[],1);
Z4                                                      = [ones(size(VarX(:)))];
Y4                                                      = VarX(:);

% Estimate coefficients
Y                                                       = [Y1; Y2; Y3; Y4];
Z                                                       = blkdiag( Z1, Z2, Z3, Z4 );
M                                                       = Z\Y;
M(5:8)                                                  = sqrt(abs(M(5:8)));
M(4)                                                    = 4*M(4);
M(3)                                                    = 3*M(3);
M(2)                                                    = M(2);
M(5)                                                    = M(5)/1.5;

end

